A quantum Monte Carlo algorithm for softcore boson systems 
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An efficient Quantum Monte Carlo algorithm for the simulation of bosonic systems on a lattice in 
a grand canonical ensemble is proposed. It is based on the mapping of bosonic models to the spin 
models in the limit of the infinite total spin quantum number. It is demonstrated, how this limit 
may be taken explicitly in the algorithm, eliminating the systematic errors. The efficiency of the 
algorithm is examined for the non-interacting lattice boson model and compared with the stochastic 
series expansion method with the heat-bath type scattering probability of the random walker. 



During the last few years there was an increasing num- 
ber of reports on strongly correlated quantum systems. 
A lot of attention has been focused on quantum phase 
transitions pj at zero temperature, which can be ob- 
served when parameters such as the particle concentra- 
tion and/or the interaction constants are varied. In order 
to observe the quantum phase transition experimentally, 
one must be able to precisely control the parameter (s), 
driving the transition, which is usually very difficult in 
real experimental situations. Therefore, only analyti- 
cal theories and numerical simulations have been able 
to provide an accurate description of the critical behav- 
ior, associated with quantum phase transitions. Quite 
recently, however, a very precise tuning of parameters 
was achieved in a system of ultra-cold atoms trapped 
in an optical lattice, formed by the intersection of laser 
beams |2j . A transition from Mott insulating phase to a 
superfluid phase was observed. It was argued that the 
system is well described by the bosonic Hubbard model 
on a d-dimensional lattice, and comparisons were made 
with numerical simulations 0,0]. This is just one exam- 
ple of an experimental realization of a strongly correlated 
quantum system, and a lot of experimental work will be 
done along these lines in the nearest future. We believe 
that it is very important in such studies to be able to pro- 
vide an accurate and simple theoretical description of the 
experimental system. Since the analytical solution of the 
models of strongly correlated systems is usually impos- 
sible, such a description may be in most cases provided 
only by the numerical simulations. 

While efficient and powerful Quantum Monte Carlo 
(QMC) algorithms exist for general quantum spin sys- 
tems, the progress in the development of the algorithms 
for the numerical simulations of bosonic systems with no 
hard core is much more modest. In the present paper 
we describe a novel QMC algorithm, allowing efficient 
simulations of the bosonic models with short-range in- 
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teractions on a lattice in the grand canonical ensemble, 
with average particle number controlled by the chemical 
potential. 

Recently an efficient QMC algorithm for the simula- 
tion of spin models with arbitrary spin quantum num- 
ber S on the lattice was proposed and implemented pj. 
It is based on coarse- graining of the conventional loop 
algorithm with split-spin representation, in which each 
spin-S 1 operator is replaced by a sum of 2S Pauli matri- 
ces. One update cycle of worldline configuration in this 
algorithm consists of a) placement of the vertices on the 
space-time lattice; b) creation of a pair of spin-raising or 
spin-lowering worms; c) propagation of one of the worms 
through the lattice with scattering on the vertices, re- 
sulting in changes of worldline configuration; d) worm 
annihilation. The algorithm for a particular model is 
thus defined by specifying a number of parameters, de- 
pending on the local worldline configuration: density of 
vertices, scattering probabilities at vertices and the prob- 
abilities for creation and annihilation of a pair of worms. 

Holstein-Primakov (HP) transformation gives a re- 
lation between the spin systems and the boson systems. 
In spin wave theories, the transformation is used for map- 
ping a spin problem into a boson problem. Here we 
do the opposite in order to derive a Monte Carlo al- 
gorithm for bosonic systems from the above-mentioned 
one for spin systems. The relation can be written as 

S+ = b\{2S - 6j& 4 ) 1/2 . ^ r = ( 2S ~ btbi) 1/2 bi> and 
Sf = ni — S, where S^',S i and Sf are spin operators 

on the site i, and b\ and bi are the boson creation and 
annihilation operators. At a first glance it appears that 
the algorithm derived from the HP transformation would 
be directly applicable only to the boson systems that have 
an artificial limitation of number of particles per site (i.e., 
it cannot exceed 25). We show that this is not the case 
in the following. 

We lift the limitation by taking the limit of large 5. 
Examining the HP transformation, we note that if n% will 
be kept finite by the chemical potential, in the large 5 
limit we can neglect higher order terms making no error, 
and keep only the lowest order in m — b\bi in the HP 



2 



transformation which leaves us with 
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2S 
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Mapping JIJ allows to rewrite the Hamiltonian of a 
bosonic model in terms of spin operators. Thus, if there 
is an algorithm for spin systems with arbitrary S, and if 
the infinite 5* limit of this algorithm exists, we can easily 
obtain an algorithm for the bosonic systems. 

To demonstrate this idea, we consider a simple model 
of non-interacting softcore bosons on a <i-dimensional hy- 
percubic lattice of linear size L with the Hamiltonian 
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where t is the (positive) hopping amplitude, fj, is the 
chemical potential and the first sum is over the pairs 
of nearest-neighbor sites. Using the mapping we can 
replace the bosonic operators in J2J with the spin opera- 
tors, leading to a model equivalent to the original bosonic 
model in the limit of infinite S: 
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Since this is an XY spin model, an efficient algorithm 
is available for any S y| Our task is, therefore, to take 
the infinite S limit of the algorithm. It turns out that 
all the parameters defining the coarse-grained algorithm 
have well-defined values in this limit as well. Below we 
describe the procedure of taking this limit and give a de- 
tailed description of the softcore boson algorithm for the 
non-interacting model. Generalizations to models with 
interactions, such as the on-site repulsive interaction and 
short ranged repulsive and/or attractive interactions, are 
straightforward. This, for instance, makes the present 
idea readily applicable to the boson Hubbard model. 

Naturally, boson occupation number must be positive, 
which leads to a restriction on the possible values of 
chemical potential: /i < — dt or > dt. To apply the 
coarse-grained algorithm we can use the values of param- 
eterSjderived for the general XXZ model in Table I of 
Ref. U Relationship between the parameters in Ref. 5 
and the parameters of our model is 
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One has to use the results of Ref. with caution, since 
they are given for the case of positive h. Therefore, in 
order to use them for the present problem, we need to 
change the sign of the field in J3J| and at the same time 
reinterpret particle numbers denoted by / and m in Ref. 
5. Namely, in the present article I denotes the number of 
holes, whereas I = 2S — I denotes the number of particles. 
Accordingly, while taking the infinite S limit with fixed 
density of particles we have to assume that I and m are 
close to 2S, whereas I and fh are of order unity. 



Probability of creation of a pair of spin-raising or 
particle-number-decreasing (PND) worms in the coarse- 
grained algorithm is 1/2S and that of a pair of lowering or 
particle- number-increasing (PNI) ones is 1/2S. By taking 
the limit S — > oo we find, that the probability to create 
a pair of PND worms is zero. Corresponding probability 
for a pair of PNI worms is then unity, indicating that our 
cycle will always start with a pair of PNI worms. That, 
however, does not mean that the number of particles will 
be constantly increasing, since the worm changes its type 
to the opposite one every time it changes direction as a 
result of scattering on a vertex. Once the traveling worm 
returns to the point of origin, it can either annihilate 
there, ending the cycle, or pass through. The probability 
of annihilation of a pair of PND worms is l/l and zero 
for the PNI ones. 

Remaining parameters, such as density of vertices and 
the vertex scattering probabilities, needed for the con- 
struction of the algorithm, can be derived by examining 
the values in Table I of Ref. H for region IV and taking 
the value of S to infinity. First of all, the vertex density 
B is given by B — h{lm + Irh + lm)/2. To list non-zero 
scattering probabilities, using the notation of Ref. |^, we 
have 
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Here we have set J' = and the superscript + or — indi- 
cates that the type of the incoming worm is PND or PNI 
respectively. Seemingly, there is a problem with density 
of vertices becoming infinite in the infinite S limit. How- 
ever, it should be noted that all non-trivial probabilities 
of the scattering events are proportional to l/B, so that 
the density of the scattering events remains finite. In 
other words, in the limit of infinite S the situation is iden- 
tical to the one that occurs when taking the continuous 
imaginary time limit in a conventional loop algorithm • 
Exploiting the analogy to the continuous imaginary time 
loop algorithm, we can easily construct a procedure for 
finding the time of the next scattering event. Namely, in- 
stead of examining each vertex, it is possible to generate 
the time of next event as a Poisson-distributcd random 
number where the average time interval or the density 
depends on the local spin configuration and the type of 
the scattering process. 

We can readily obtain the density of such events by 
multiplying the scattering probabilities ©, © an d 
by B and take the infinite S limit. Since Eq. J7J) yields 
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zero, we have two non-zero scattering densities for inter- 
vals that have no kinks in it: 
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and 
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For the scattering probability at kinks, only the scat- 
tering probability (JSJ) will remain non-zero , since the 
probability © vanishes in the infinite S limit. 

In order to describe the algorithm in detail, we intro- 
duce a concept of a constant environment interval (CEI) 
on which the moving worm resides. A CEI is defined as 
an interval ahead of the worm in which the environment 
of the worm does not change in the imaginary time di- 
rection, i.e., the worldline state changes neither on the 
current site nor on any of the neighboring sites. This 
interval is bounded by one of three events, closest to the 
worm: a) a kink on the current site, b) a kink on one of 
the neighboring sites, or c) the point of origin, where the 
other worm waits for the moving worm. 

Worldline configuration update cycle for the non- 
interacting model may be summarized as follows: 

1) Choose an arbitrary space-time point to place a pair 
of worms, one of which will move, producing the changes 
in the configuration, and another one will mark the point 
of origin. Always start with a PNI (spin-lowering) worm. 
Choose the arbitrary direction (up or down) for the 
worm's initial movement. 

2) Determine the CEI. 

3) For each type of scattering and for each nearest neigh- 
bor site, which is a candidate for the final scattering des- 
tination, generate the time of the next possible scattering 
event stochastically according to the Poisson distribution 
with the densities (JTDJ and lfTT)l . 

4) If the advancement of the worm by the smallest of 
these times does not take the worm out of CEI, imple- 
ment the corresponding scattering event. In case of a 
back-scattering event, change the type of the worm to 
the opposite one. Go back to 2. 

4') If the advancement of the worm get the worm out of 
CEI, advance the worm to the end of the CEI. 

5) If the end point of the CEI is not a kink or the original 
starting point, go back to 2. 

5') If the end point of the CEI is a kink, attempt to scat- 
ter on it according to JSJ. If the scattering fails or not 
applicable, let the worm skip the kink and go on. Go 
back to 2. 

5") If the end point of the CEI is the original starting 
point, stochastically determine, whether they will anni- 
hilate with the probability l/l where I is the particle 
number on the CEI. If it annihilates, the update cycle 
is terminated. Otherwise go back to 2. 
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FIG. 1: The compressibility plotted against the average oc- 
cupation number for three-dimensional free lattice boson sys- 
tem at ksT/t — 2.0. The lines are the exact analytical values, 
while the symbols are the results of QMC simulation. Uncer- 
tainty in the occupation number is in all cases smaller than 
the width of the symbols. Data for L — 16 and (n) > 0.4 
is computed with the standard number of cycles (solid dia- 
monds), while the data for (n) < 0.4 was obtained with a 100 
times longer simulation (open diamonds). Error bars are in 
all cases smaller than the size of the symbols. 



After a number of full update cycles (resulting in worm 
annihilation) the observables are measured. 

To test the validity and evaluate the efficiency of the 
algorithm we have performed a number of tests, compar- 
ing the results of QMC simulation of the non-interacting 
boson model in three dimensions to the exact results. Al- 
though the model does not have any interaction terms, it 
is non-trivial enough to provide us with excellent grounds 
for testing because it displays Bose-Einstein condensation 
(BEC) and the observables may be calculated analyti- 
cally. 

We have performed simulations at k^T — 2t at 10 dif- 
ferent values of the chemical potential, chosen so that 
the resulting average occupation number would be n = 
(N) /V = 0.1, 0.2, . . . , 1.0. Three different system sizes 
are considered: L = 4, 8, 16. If not stated otherwise, 
for each value of system size and chemical potential we 
have performed 50,000 cycles for equilibration, and an- 
other 50,000 cycles for measurement. The 50,000 mea- 
surements cycles were divided into 10 bins of the equal 
length for estimating the statistical error. 

In all cases we investigated, including cases close to 
criticality and ones deep inside the superfluid phase, 
we found an excellent agreement between the numeri- 
cal QMC data and the exact analytical results. As an 
example, Fig. ^ shows the dependence of the compress- 
ibility K = {dn/dfi) T = (fc B r) _1 ((A r2 ) - (N) 2 )/L d as 
a function of n at a fixed temperature k^T = 2t. The 
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FIG. 2: The superfluid density plotted against the average 
occupation number for three-dimensional free lattice boson 
system. Standard simulation parameters were used for all 
data points. The lines are the exact analytical values while 
the symbols are the results of QMC simulation. 

observation of the divergent behavior of compressibility 
at BEC transition is well within the reach of numerical 
simulation. For low values of n and L — 16 we had to in- 
crease the number of cycles 100 times in order to obtain 
a good statistics because the typical lifetime of a worm 
becomes too short in this case. Even after this increase, 
the CPU time spent for this case is smaller than that for 
the superfluid cases. 

In Fig. [21 we plot the superfluid density ps against 
the average occupation number, ps is defined [1] as 
p s = L- d k B T[d 2 F{6)/de 2 ] e=0 where F{6) is the free 
energy of a system twisted by the angle 6 per lattice spac- 
ing. In QMC simulation this quantity can be measured 
by p s = L 2 - d k B T(W 2 } where W x is the sum of winding 
numbers of all world lines in the x-direction. The possi- 
bility of measuring the winding number fluctuation is one 
of the advantages of the present approach, compared to 
the algorithms, such as the one used in Ref. 0, that works 
in the fixed winding number ensemble. In Fig. [21 we can 
again see that the onset of the condensation is captured 
by the QMC simulation with the present algorithm. 

We have also compared the performance of the new al- 
gorithm with the directed loop algorithm [9} , one of the 
best QMC algorithms currently available for the simula- 
tion of softcore boson systems. The directed loop algo- 
rithm is quite general and powerful method, applicable, 



in principle, to any quantum system. However, it is up 
to the user to find a set of scattering probabilities, opti- 
mizing the efficiency of the algorithm for a given model. 
Due to a huge freedom in the choice of algorithm pa- 
rameters, in most cases this is a highly nontrivial task. 
In addition, to apply the directed loop algorithm to the 
softcore boson systems, one has to set an artificial up- 
per bound for the site occupation number. This upper 
bound must be taken large enough to make the simu- 
lation free from the systematic error. This could be a 
serious disadvantage especially for bosonic system with 
a random chemical potential where a large number of 
particles may reside on the same site. The present algo- 
rithm is free from these disadvantages. For comparison 
purposes we have used the directed loop algorithm with 
a set of simple heat-bath scattering probabilities and the 
site occupation number was limited by rii < 20. 

While it is hard to make a quantitative comparison of 
two essentially different algorithms, we have established 
that our proposed algorithm in all considered cases per- 
forms better, i.e. in all cases smaller computational time 
was needed to reach the same error bars. While the com- 
parison of computational times does not inequivocally 
prove that our algorithm is more efficient, we note that 
in the vicinity of the critical point and in the superfluid 
phase we have failed to obtain the reliable estimates for 
the observables with the directed loop algorithm within 
reasonable computational time. 

In summary, we have described a construction of an 
efficient QMC algorithm for the simulation of softcore 
boson models on the lattice, based on the coarse-grained 
algorithms for the spin models. By establishing the re- 
lationship between the boson and spin operators in the 
infinite S (total spin quantum number) limit, we have 
mapped the model of non-interacting bosons on a lattice 
to a spin XY-mode\ in a magnetic field. We have demon- 
strated, that the limit of infinite S may be taken directly 
in the algorithm, leading to improved performance and 
absence of systematic errors. The resulting algorithm 
was found to perform better than existing algorithms. 
The result of applications of the present algorithm to 
other models, such as the Bose Hubbard model, will be 
reported elsewhere [Tof . 
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